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We investigate the nature of the phase transition occurring in a planar XY-model spin system 
with dipole-dipole interactions. It is demonstrated that a Berezinskii-Kosterlitz-Thouless (BKT) 
type of phase transition always occurs at a finite temperature, separating the ordered (ferro) and the 
disordered (para) phases. The low-temperature phase corresponds to an ordered state with thermal 
fluctuations, composed of a "gas" of bound vortex-antivortex pairs, which would, when considered 
isolated, be characterized by a constant vortex-antivortex attraction force which is due to the dipolar 
interaction term in the Hamiltonian. Using a topological charge model, we show that small bound 
vortex pairs are easily polarized, and screen the vortex-antivortex interaction in sufficiently large 
pairs. Screening changes the linear attraction potential of vortices to a logarithmic one, and leads to 
the familiar pair dissociation mechanism of the BKT type. The topological charge model is confirmed 
by numerical calculations, in which we demonstrate that the transition temperature slightly increases 
when compared with the standard BKT result for short-range interactions. 



The standard paradigm of phase transitions in a pla- 
nar system of electrically neutral particles representing 
an effective XY-model spin dictates that long-range or- 
der does not exist at any finite temperature [1, 2J. On 
the other hand, below a critical temperature Tc, spin- 
spin correlations in two spatial dimensions decay in a 
power-law fashion. As a consequence, short-range co- 
operative phenomena, such as superfluidity, can exist 
at these low temperatures. Below the critical point, 
thermal excitations predominantly occur in the form of 
vortex-antivortex pairs (VAPs). Due to the attraction 
between vortex and antivortex, the VAPs remain bound 
at T < Tc, and the gas thus remains superfiuid, while 
the dissociation of VAPs at T > Tc leads to an exponen- 
tial decay of correlations. This standard scenario, corre- 
sponding to the celebrated BKT phase transition |3H5], 
however conventionally only applies when the particles 
interact by short-range (contact) interactions. 

Spin ordering phenomena in planar systems with long- 
range forces mediated by the dipole-dipole interaction 
(DDI) between particles carrying static dipole moments 
/i, have been extensively studied in the past, cf., e.g., [6- 

. It was found that the DDI leads to a stabilization of 
the long-range ferromagnetic order against thermal fluc- 
tuations, and the ground state of the spin system may 
thus be spontaneously polarized |6HT2] , or acquire various 
structures |14fl^ . So far most efforts aimed at under- 
standing the phase transition (PT) were a combination of 
renormalization group (RG) arguments with phenomeno- 
logical approaches [7^151 , or Monte Carlo and Molecular 
Dynamics simulations [19-24_. However, it is fair to say 
that a complete picture of the nature of the PT in the 
dipolar XY model so far has been elusive. There is, in 



particular, no agreement on both the physical mecha- 
nism of the PT in the dipolar XY model and its critical 
temperature. 

For contact interactions, vortex and antivortex in a 
VAP interact by a logarithmic interaction potential, 
which in superfluids is due to kinetic energy of the flow, 
and the attraction force decays with the inverse distance. 
This logarithmic interaction is the primary requirement 
for the occurrence of the BKT transition. On the other 
hand, in the presence of a dipolar interaction term in the 
Hamiltonian, a spatially constant attraction force Kq be- 
tween vortex and antivortex occurs jH HH HI] (also see 
below), which bears an obvious potential importance for 
the PT, which has been overlooked in most previous in- 
vestigations of the dipolar XY model. We note that the 
confinement in isolated VAPs by the string tension is 
one of the rare instances outside the realm of Quantum 
Chromodynamics (QCD), in which linear interaction po- 
tentials between a particle and its antiparticle, in QCD 
between quark and antiquark, occur. 

The existence of vortex conflnement led to the state- 
ment HH [12 that a novel PT distinct from BKT takes 
place in the two-dimensional (2D) dipolar gas. In the fol- 
lowing, we critically examine this conclusion. We demon- 
strate that the sole effect of the DDI is that a VAP dis- 
sociation transition of the BKT type occurs at slightly 
higher temperatures. To this end, we use an analyti- 
cal model, applicable sufficiently close to the transition 
point, backed up by numerical simulations for the full 
range of temperatures. We show that the apparently in- 
consistent pictures of conflnement of isolated VAPs and 
occurrence of a PT driven by the familiar BKT mech- 
anism of pair dissociation can be made fully consistent 
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with each other if one correctly accounts for shielding of 
the bare VAP tension at finite temperatures. The shield- 
ing effect of the linear attraction potential within a large 
pair is mediated by the large number of small VAPs in 
which it is immersed. The vortex-antivortex interaction 
only remains linear at small distances R between the vor- 
tices in a pair, R < tq, where the length parameter tq is 
defined in Eq. ([t]) below; it however becomes logarithmic 
at sufficiently large VAP size, R > tq. 

The polarization (effective spin) states of the dipo- 
lar XY model, which correspond to a VAP gas, are 
described by the continuous two-component vector-field 
s(r) — scosBex + ssinOey, representing the hydrody- 
namic, coarse-grained average of the spin, where 9 is the 
azimuthal angle. In dimensionless form, the hydrody- 
namic free energy functional is Gs [s{r)] — Go + Gdd- 
The short-ranged term Go corresponds to a perturbation 
of the uniformly polarized ground state (a, b = x,y) [27\ : 



(1) 



where df is the surface element, while the dipolar inter- 
action energy functional is given by [SI [71 [TTJ [H] 



Gm = A dfdf 



,Pp ir)pp {r') 



(2) 



Here pp (r) — —V • s (r) is the density of polarization 
charges, and A oc /i^ a dimensionless interaction constant 
(representing the DDI coupling constant). The appropri- 
ate units of energy and length depend on the specific sys- 
tem considered; as its value does not affect our argument, 
we have also put the "contact interaction" coefficient of 
the (s^ — 1)^ term equal to unity. We assume that the gas 
of dipoles is stable against collapse due to the long-range 
attractive part of the DDI (head to tail instability) at 
any value of A because of a sufficiently large short-range 
repulsive core of the true interaction potential. Exam- 
ples for possible applications of the present model span a 
broad range, including the commonly studied thin mag- 
netic |14j and ferroelectric films, quantum gases [28-r 
|3T| . but also, in a biological context, water molecules on 
hydrated proteins |32J. 

Focusing our attention on the regime T sa Tc, we have, 
typically, that the VAP size i? > 1. For these VAPs, 
large compared to a microscopic (e.g. coherence length) 
scale, using that V • s 1/R and hence Go ^ log-R, 
Gdd ^ R- Therefore, the energy of a VAP with R ^ 1 
increases linear with their size, and vortex and antivortex 
attract each other with a constant force [H [TTJ [12], in 
sharp distinction to contact interactions, where the force 
decreases linearly as the VAP size increases. We then 
restate the total free energy of the system in the form 



(i, j = 1 . . . 2N, where N is the number of VAPs) [TTIIT^ 

The function u [r) w 2it log (1 + ar) + Kqt, and qi = ±1 
are the topological charges of the vortices |1! IS] located at 



\ri — rj\. The quantity Kq = Kq (A) 



is the "VAP tension" coefficient, depending on the cou- 
pling strength A. The logarithmic contribution in u (r) 
arises from Go in Eq. ([ij and describes the interaction 
of vortices in the limit of A — >■ |3H5]. The parameter 
a describes interaction between vortices at r — >■ and 
when A = 0; it is found from numerical calculations of 
the VAP energy (see below). At r — > vortices of a 
given pair annihilate, corresponding to the requirement 
u (0) — 0. The "potential" $ (r) in ^ is introduced in 
analogy with electrostatics: 



$ (r) = ^ QjF (r - rj) = J F (R) p (r') d^r'. 



(4) 



Here, R ^ r - r', F{R) ^ -27r In (1 + ai?) - KoR, 
and the vortex topological charge density p (r) — 
J2j QjS^^^ {'"' ~ '"'])■ The topological charge inside a 2D 
contour G is defined as usual to be Q = (27r)~^ dO. 

Eq. ^ is the central assumption leading to the phe- 
nomenological description of our model system. We have 
set up the corresponding Langevin Dynamics simulations 
and ran extensive checks to establish the (nontrivial) va- 
lidity of Eq. ([3]) (see supplementary material) . Numer- 
ical calculations also let us establish the functional de- 
pendence of ivro(A). We found that at A > Kcr ~ 0.4, 
corresponding to (Acjj) w 3.1, a physical instabil- 
ity of the single VAP configuration arises: A new small 
VAP is spontaneously created in the center of a large 
VAP. With further increase of A, the "parent" vortices 
are immersed into a cloud of small polarized VAPs that 
decrease the effective topological charge of the parent 
vortices. According to Q, the topological charge de- 
crease is equivalent to it'o increasing with A, finally sat- 
urating at Kq (A) « Kq (Aci?) « 3.1 (see the inset of 
Fig. [2]). The spontaneous creation of pairs follows also 
from the expression ([3]) : At A > A^^ the energy of the 
large (i? ^ 1) pair decreases with the emergence of a 
new small pair in the center. This occurs at A = ^cRi if 
we choose a — Kq [Kcb) /2tt ~ Kq (oo)/27r. Typical nu- 
merically calculated (see supplementary material) vector 
field s (r) configurations and total energy distributions 
at r = and close to Tc are depicted in Fig.[T] illustrat- 
ing the screening effect for the linear attraction by the 
small VAPs close to the PT. 

Once the additive expression for the energy ([3| is 
verified by the numerics, we can use it to elucidate 
the nature of the PT and to derive an analytical ex- 
pression for the transition temperature T^. By its 
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Figure 1: Left: Polarization field s(r) and total energy dis- 
tribution at T = for a VAP located at (0,±30). Brighter 
regions characterize larger energy density, where vortices or 
VAPs are located. The angle between the axis of the large 
VAP and the asymptotic uniform polarization is chosen to be 
45°. Right: Configuration at temperatures close to Tc, where 
a large number of thermally excited small VAPs strongly al- 
ters the power law of attraction in the given VAP at (0, ±30). 



definition, the potential $ (r) satisfies a Poisson type 
equation, (r) = p{i"). Here the Hnear operator 
Lr is defined such that LrF {r — r') = J*^^^ (r — r'), 
which in Fourier representation reads Lk = 1/Fk = 
{2TTKok-^ + 47r2a/ [k^{k + a)] y\ Let us follow the 
electrostatic analogy further. Since the energy of a charge 
(vortex) q placed in the external potential $ is J7 = 5$, 
the force, acting on the charge (vortex) is = qE, where 
the quasi-electric field vector E — — V$. The energy of 
the pair in this field is [/ = (7+$ (r+)+(7_$ (r_) ss —d-E. 
Here the topological dipole moment of a pair is intro- 
duced according to d = qjVj = q+rj^+q^r- = r, 
where r = r_|_ — r_ . The density of the topological polar- 
ization charges is pTp — —V • P, where P (r) — np (d) 
is polarization vector of the VAP gas, and np is surface 
density of VAPs; (...) denotes statistical averaging. Af- 
ter averaging the polarization topological charge inside a 
contour C, Qtp — / dfpxp, becomes a fractional number 
in general, while before averaging it must be integer. 

The equation for the potential of a point charge Q 
placed at the origin is <i>Q (r) is Lr^g (r) = QS'^^^ (r) -|- 
Ptp- In the Weak Field Approximation (WFA), that is 
keeping only the linear terms in V<i> of Eq. ([s]), we obtain 
a general formula for P (r): 



P[r) = - j dfi; (|r - r'|) V$ (r') , 



(5) 



Hence {pTp)k = —k'^-ip(k), where ip {k) ~ 
J dfexp{—ik-r)ip{r). Since i/j (r) decays fast enough 
with r, then % = ?/> (fe = 0) = f dftp (r) is a finite number, 
which leads to {pTp)k ^ ^xk^^Qk at fc 0. Note that 



for the slowly varying field P (r) « = — xV<I> (r). 

Therefore x is the susceptibility of the VAP gas. Putting 
everything together, we have for the potential 



^Q{r) = Q 



<Pk exp {ik ■ r) 



xk"^ 



We define the distance scale 



1 



2^KoX 



(6) 



(7) 



which involves the VAP tension Kq and the suscepti- 
bility X of the VAP gas. At large distances, r ^ tq, 
the dominant contribution to the integral comes from 
small values of fc, for which the term in the denom- 
inator is negligible. Therefore, since (Inr)^ — — 27r/fc^, 
the potential of the "charge" at large distances is loga- 
rithmic: $Q(r) = — ((5/27rx) In (r/Ci) , r ^ rp, where 
the constant Ci ^ tq. In the opposite limit, the po- 
tential is linear: $q(t) « —QKqt, r ^ tq. We 
propose a simple interpolating expression: '&q(t') ~ 
— (Q/27rx) In (1 + r/rg), so that the energy of a suffi- 
ciently large pair of size R is given by 



U {R) ^ -^Q=i{R) 



1 

27x 



In 1 



R 

To 



(8) 



The distance scale rg in Eq. (|7| represents an analogue of 
the Debye screening radius for 2D interactions of topolog- 
ical charges. Indeed, the polarization topological charge 
density of VAP gas next to the single topological charge 
Q equals: pxp — —V • P ~ —xQKo/r. We conclude 
that the total charge inside a circle of radius R is given 
hy Qt = Q + Qtp « Q (1 — 2'itKqxR)- From this quali- 
tative consideration we thus come to an important con- 
clusion: the charge is essentially shielded, Qt « 0, which 
occurs at a distance scale R ^ vq. The polarization of 
the VAP gas inhibits the linear attraction within large 
VAPs, which would prevail with shielding not taken into 
account. Hence we conclude that the PT associated with 
the dissociation of pairs is qualitatively very similar to 
the BKT transition in a system with A — 0. 

The standard calculation of the transition temperature 
for a gas of polarizable VAPs with interaction ([s]) gives 
the following implicit equation for the transition temper- 
ature in terms of the susceptibility [3H5], 



T, 



1 

4ttx' 



(9) 



At the transition temperature, T = Tc, VAPs begin to 
dissociate. This implies that at Tc — T <C Tc only a 
small fraction of the pairs is large and close to dissocia- 
tion. For this reason it is possible to neglect the interac- 
tions between the largest VAPs and calculate the energy 
of a single large pair approaching its dissociation limit, 
which is permeated by a cloud of comparatively small 
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bound VAPs. Hereinafter we will subdivide VAPs into 
two classes: small pairs (SP) and large, close to dissocia- 
tion, pairs (LP). The shielding effect arises due to the po- 
larization cloud provided by SPs, influencing the poten- 
tial energy of a LP. Let us calculate first the polarizability 
ap of a single SP. The energy of the SP in an external 
field E equals V (r) = u{r) ~ r • E. The average dipole 
moment of the SP and the susceptibility of the VAP gas 
within the framework of the WFA are given by the rela- 
tions (d) = J d/r-exp i~V/T) / J d/exp (-V/T) « apE, 
X = apnp. Here, ap — I3/ {2TIi) is the SP polarizabil- 
ity we are looking for, /„ = 2n drr'^e^''^ (1 + ar) ^ , 
^ = 27r/r, 7 = Kq/T. The potential energy of a SP in 
the field of the charge Q, (d) V$ « (d) V^vip), approx- 
imately does not depend on the position of the pair, p, 
so that np can be considered as a constant. 

Using Eq. (|8| , the average interaction between 
the vortices at finite temperature is {U (i?)) = 
(2^x)"'ln(l + i?i3p/ro) - Zp^ dRU {R) g {R) - 
re(l-4r)/[(-T)(l-2r)], where zp = J^dRgiR), 
g {R) = Rexp[~U (R) /T], and the relative tempera- 
ture T = {T- Tc) /Tc. Therefore, next to the PT, 
|r| ^ 1, the dimension of a dissociated pair is large, 
Rdp ^ tq exp (1/ |t|), and its topological "electric" field 
is small, E = —VU (i?) oc exp (—1/ |t|). A condition for 
the applicability of the WFA is therefore the existence of 
a small parameter exp (—1/ |t|), representing the ratio of 
the typical topological "electric" field, ^ \VU {R)\ to the 
field inside the SP, - \Vu{R)\. Therefore, the WFA is 
applicable close to the transition temperature. 

The surface density of VAPs, rip, at Tc can be calcu- 
lated as follows. At T « Tc the VAPs only start to disso- 
ciate and the fraction of LPs is small. The typical pair is 
small, and the interaction w u (r) between its vortices is 
still unscreened. The partition function of N pairs on a 
surface with area S is given by z = z^ /N\ sa {zie/N)^ , 
with zi = f d^ pd^r exp {—u/T). Minimization of the free 
energy Fp (N) = — Tlogz gives N — zi, np — N/S = /i. 
Using this result and the relation Tc = (Attx)^^, there 
follows 2ttIj ~ 1, leading to an implicit equation for the 
critical temperature 
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Figure 2: Transition temperature Tc versus DDI coupling 
strength A computed using Eq. ( |10[ ) (dashed line), and nu- 
merical result for Tc (diamonds). In the inset the saturation 
effect for the VAP tension coefficient Kq with increasing DDI 
strength A is shown. The functional dependence Kq = Ko{A) 
forms the input for the evaluation of Eq. ( 10 1 for the critical 
temperature Tc- 



47r" 



dr 



r3 exp i-Kpr/Tc) 
(1 + ar)'"/^= 



(10) 



we derive Tc (00) « 1.3. Note that the RG arguments 
of jTTl [12] gave the much larger prediction Tc{oo) — 2tt. 
This in turn implies that it is rather difficult to account 
for all essential Feynman graphs in order to adequately 
describe the shielding effect in a RG calculation, which 
applies equally well to the case A = 0. 

In conclusion, we have demonstrated both numerically 
and analytically using an analogy to plasma physics, that 
VAPs in the 2D dipolar XY-model dissociate at a critical 
temperature in a manner familiar from the BKT transi- 
tion. This is due to the linear confinement potential in an 
isolated large VAP being shielded by a gas of small VAPs, 
in which the large pair becomes immersed around the 
transition point. Therefore, the logarithmic attraction 
between vortices in large pairs is restored. By obtaining 
a physically transparent scenario, we have therefore pro- 
vided an unambiguous proof that the BKT mechanism 
is applicable to a much broader class of systems than 
hitherto established. 
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Next we compare the semianalytical solution for Tc 
from ( p!o| with the results of the numerical calculation, 
cf. Fig]2j As described in detail in the supplementary ma- 
terial, we used Langevin Dynamics and Binder's method 
to calculate Tc in a series of simulations of increasingly 
larger realizations of the model system. In the non- 
interacting case, Kq = (or A = 0), the equation (10) 
yields Tc ~ 0.7. From our numerical calculations, we 
found that Tc = 0.85. At large A the transition tempera- 
ture tends to the constant value Tc (00) « 1.4. From (10) 
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by: 

^ = -7^(Go + G,,) + Ca, 
at asa 

where 7 is a constant that determines the time scale of 
relaxation. The stochastic thermal noise terms satisfy 
iUt)) = and iUtW)) = '^T^S^pSit - t'). In the 
LD simulations we use a second-order Runge-Kutta al- 
gorithm. The equations of motion above are integrated 
numerically with the sufficiently small discrete time step 
At = 0.005. To compute the DDI term ^ in an efficient 
0{Ng\nNg) way we used a NumPy FFTW realization 
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Figure 3: Free energy of a uniformly polarized 2D spin system 
with a small VAP in the center of a large VAP (a) and for a 
slightly shifted small VAP (b) versus the angle (in units vr/12) 
between the axis of the small pair and the x axis. Red (N) 
and blue (A) curves correspond to the numerical calculation 
and the analytical approximation ([3|, respectively; A = 0.2. 



Supplementary material 

We studied the thermodynamics of our model us- 
ing Langevin Dynamics (LD) [SlJ . Using a discretized 
representation of the Hamiltonian at every grid point 
a = 1,2, Ng, we performed a fixed temperature run 
of sufficient length to get reliable averages. The calcula- 
tions were performed using periodic boundary conditions 
on square lattices L x L with the number of independent 
nodes Ng = L^. The LD dynamical equations are given 



The software was first used to check the assumptions 
leading to Eq. (|3| . First we consider the case of zero tem- 
perature, r = 0; and checked the additivity rule ([3| nu- 
merically by using the imaginary time relaxation method, 
which is capable to find local minima of the free energy 
functional on a configuration space of 2D vectors Sa taken 
in nodes of a dense square lattice. Following |S3| . the ini- 
tial approximation was specified by the complex-valued 
skyrmion type expression sq{z) ~ 2W/ ^1 + , with 

W = Y[^=i ~ ^-jf^ ; ^'^d z = a; -f «y, Zj = xj + iyj. 
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merically the total energy within a small error, less than 
~ 5% (examples of typical distributions are presented in 
Fig.|3|. The string tension constant Kq was calculated 
by analyzing the forces, acting on vortices in a VAP by 
calculating averages of the derivatives of the pinning po- 
tential with respect to the positions of the vortices. The 
accuracy we were able to achieve is limited by the per- 
turbations introduced by the pinning potential, and is 
sufficient to prove the reliability of the additivity rule 
Eq. (|3|, also cf.Fig.jsj 



Figure 4: Explanation of Binder's approach to obtain the crit- 
ical temperature S4J. Binder's parameter Ul was calculated 
for A — 0.2; the intersection point for different lattice sizes L 
yields (indicated by the bold arrow). 



For a gas of pairs the total charge vanishes, J2j Qj = 0- 
Therefore, — )■ 1 at z — )■ oo, which agrees with the 
physical boundary condition s — (1, 0) at r — )■ oo. At the 
first stage of imaginary time propagation the cores of the 
VAPs begin very slowly to approach which, finally, leads 
to their annihilation. Due to this mutual attraction, the 
VAP therefore is not a real local minimum. Any initial 
configuration inevitably transfers at T = to the abso- 
lute ground state s — >■ (1,0). To stabilize VAPs, "pin- 
ning" of vortices was used by adding to the free energy 
a term Gpi„ = G,, G, = f dfVj (r) [s (r) - So (r)]^ 



1. We in- 



where Vj (r) = Vo cxp — (r — rj ) /a 
vestigated multiple configurations with different numbers 
of VAPs using the expression ([s]), and reproduced nu- 



To explore the critical behavior of the model depend- 
ing on the DDI coupling constant A numerically we 
use Binder's method |S4| and calculate the parameter 
Ul ^ I - (s^)/3/(s2)^ ("Binder's cumulant") versus 

temperature. Here s = S/Ng, S = X]^=i *a ^^e sta- 
tistical averaging is done in a manner equivalent to the 
average over polarization configurations obtained with 
LD. The intersection point of Binder's cumulants for dif- 
ferent values of the system size L gives T^, as shown in 
Fig.g 
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